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ABSTRACT 


The  use  of  adaptive  control  algorithms  was  studied  for  microprocessor  driven 
direct  digital  control  of  elementary  heating  and  cooling  subsystems.  An 
algorithm  was  designed  for  digital  regulation  of  a linear,  time-invariant 
first-order  system  with  a system  dead  time.  A recursive  least  squares  algor- 
ithm was  used  to  estimate,  on-line,  the  parameters  of  the  time- invariant  linear 
system.  The  parameter  estimates  were  then  used  to  calculate  the  feedback  gains 
of  a Proportional  plus  Integral  (PI)  controller. 

Through  computer  simulations,  the  adaptive-parameter  Pi-controller  was  compared 
with  a constant-parameter  Pi-controller . On  the  basis  of  favorable  simulation 
results,  the  adaptive  algorithm  was  implemented  for  direct  digital  control  of 
an  air  handling  unit  in  a laboratory  building  at  the  National  Bureau  of  Standards, 
Gaithersburg,  Maryland.  The  convergence  of  the  parameter  estimates  and  the 
step  response  proved  to  be  satisfactory  provided  the  system  was  operating  in  a 
linear  or  weakly  non-linear  region,  and  was  in  steady  or  quasi-steady  state. 

By  selecting  a proper  scale  factor,  improved  performance  may  be  obtained  when 
system  characteristics  vary. 

Key  words:  adaptive  control;  air  handling  unit;  direct  digital  control; 

energy  management  and  control  systems;  HVAC  system  control; 
parameter  estimator;  Pi-controller ; recursive  least  squares 
algorithm;  self-tuning  control  algorithm 


* Ur.  David's  contribution  to  this  research  effort  was  carried  out  while  he  was 
employed  by  NBS  from  April  1980  to  May  1981. 
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1.  INTRODUCTION 


A building  system  requires  an  integrated  system  of  controllers  to  provide  the 
best  possible  level  of  comfort  with  a minimum  level  of  energy  consumption. 

One  well-known  type  of  controller  is  the  Proportional  plus  Integral  plus 
Derivative  (PID)  controller.  Conventionally,  PID  controllers  are  tuned  man- 
ually to  achieve  accurate  performance.  Unfortunately,  this  manual  process  is 
not  a simple  task  because  of  many  factors  to  be  accounted  for,  and  moreover, 
a controller  properly  adjusted  for  one  season  is  not  likely  to  perform  properly 
at  a different  time  of  year.  This  problem  suggests  the  use  of  modern 
estimation  and  control  algorithms  for  performing  the  tuning  process  in  real-time. 

The  benefits  of  using  microprocessors  for  tuning  controllers  have  been 
recognized  in  many  industries  [1-4].  The  parameters  of  a controller  can 
be  adaptively  updated  based  on  estimates  of  the  parameters  of  the  system 
to  be  controlled  [5-7]  using  the  methods  of  modern  control  theory.  There 
are  two  recognized  classes  of  adaptive  control  methods — the  model  reference 
method  [8]  and  the  self-tuning  regulator  [9]. 

Kurz,  Isermann,  and  Schumann  [10]  have  compared  various  parameter-adaptive 
control  algorithms.  They  indicate  that  parameter-adaptive  control  can  be  used 
for  tuning  of  digital  control  algorithms,  for  digital  adaptive  control  of  slowly 
time-varying  processes,  and  for  digital  adaptive  control  of  weakly  non-linear 
and  partially  unstable  processes. 


For  our  purposes,  the  self-tuning  regulator  approach  seemed  most  appropriate. 

The  study  was  begun  by  deriving  a model  of  a simple  heating,  ventilation 
and  air  conditioning  (HVAC)  system  and  using  this  model  to  simulate  the 
performance  of  a control  system  consisting  of  a recursive  least  squares 
algorithm  for  parameter  identification  and  an  adaptive  Pi-algorithm. 

Based  upon  the  simulation  results,  an  algorithm  was  developed  for  direct 
digital  control.  In  laboratory  experiments,  the  algorithm  was  employed  in 
microcomputer-based  direct  digital  control  of  an  air  handling  unit. 

2.  MATHEMATICAL  MODEL 

For  simulation  purposes,  we  modeled  an  air  handling  unit  as  shown  in  figure  1. 
Outdoor  air  enters  the  system  through  the  opening  of  dampers  and  is  mixed  with 
return  air.  The  mixed  air  passes  through  an  air  filter,  a heating  coil,  and 
a cooling  coil.  A fan  supplies  the  conditioned  air  to  the  demand  zone  [11]. 

In  this  simple  system,  the  mass  flow  rate  of  hot  steam  is  controlled  by  a 
proportional  valve,  which  is  governed  by  a controller.  The  controller  operates 
to  minimize  the  difference  between  the  temperature  measured  by  a sensor  in  the 
supply  duct  and  the  set  point  value.  The  flow  rate  of  chilled  water  is  assumed 
to  remain  constant.  For  simplicity,  feedback  by  a room  thermostat  is  not 
considered.  The  problem  addressed  is  to  find  algorithms  for  tuning  conventional 
digital  controllers  by  providing  appropriate  system  parameter  values.  Thus,  we 
replace  the  conventional  controller  with  an  adaptive  digital  controller 
implemented  via  a microcomputer. 

Following  the  self-tuning  regulator  methodology,  we  depict  the  interrelationships 
among  the  elements  of  the  control  system  in  figure  2.  The  control  system  consists 
of  a process,  a valve,  a controller,  and  an  estimator.  The  process  produces  the 
output  y using  the  change  of  valve  output  Auv  as  an  input.  The  process  output 
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then  combines  with  outdoor  and  return  air  disturbances  w resulting  in  y.  The 
difference  between  the  resultant  output  and  the  set  point  value  (reference  value) 
r yields  an  error  e.  When  parameters  of  the  controller  are  fixed,  the  controller 
output  uc  depends  on  the  error  signal.  But  if  the  controller  is  adaptive,  it 
depends  not  only  on  the  error  e but  also  on  the  new  estimates  of  0 of  the  system 
parameters  as  obtained  by  the  estimator.  The  estimator  requires  the  current 
input  value  uc  and  output  value  y as  its  inputs  as  well  as  prior  values  of  uc 
and  y in  order  to  update  the  parameters.  Details  are  given  below. 


2 . 1 Process  Simulation 

We  consider  the  heating  and  cooling  coils  as  the  process,  assuming  that  the 
chilled  water  going  through  the  cooling  coil  has  constant  temperature 
and  flow  rate.  Under  this  assumption,  we  simulate  the  heating  coil  as  the 
process,  treating  it  as  a continuous  linear  device,  with  Laplace  transform 


°H(S> 


S + (X 


(1) 


where  K,  d,  and  a are  the  gain,  dead  time  (transport  lag),  and  inverse  time 
constant  of  the  system,  respectively.  The  time- invariant  values  of  K,  d,  and 
a are  determined  prior  to  simulation. 


Assuming  that  the  control  signal  is  changed  discretely  every  T'  = nT  time  units, 
where  T is  the  simulation  sampling  interval,  the  system  model  must  incorporate 
a zero  order  hold  in  cascade  with  the  coil.  The  transfer  function  of  the 


sampler  is 


-nTs 


From  equations  (1)  and  (2),  the  transfer  function  of  the  process  with  the 
and  hold  yields  in  the  frequency  domain 


H(s)  = G ( s)  G (s) 

n S 


or 


H(s)  = 


"iiTSv  /„  -dsN 
(1  - e )(Ke  ) 

s(s  + a) 


Equation  (3)  can  be  written  as 


H(s)  - [e^2  - a-(nT  + d)s](K)(i)(- 


1 + 


l' 

a • — 
s 


Using  the  block  diagram  convention  [12],  we  depict  equation  (4)  in  figure 
where  x^  and  x^  are  state  variables.  From  this  block  diagram,  we  obtain  < 
set  of  state  equations  based  on  the  dynamic  linear  network  theory.  These 
equations  in  the  time  domain  are 


y - x1 

X]_  = x2  - ax]L 
*2  = Ku 

where 


u ( t ) = u (t  -d)  -u  (t  -d  -nT) 

V V 


In  a matrix  form,  we  obtain 


sample 


(3) 


(4) 
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x = A x + B u 


(5) 


y = C x 


where 


x = fx1,  x2]  , x = l x^ , x2J 


A = 


(6) 


C = [1  0] 


In  equations  (5)  and  (6),  u and  y are  scalar  quantities  since  the  system  is  a 
single  input  single  output  system.  Given  a dynamic  linear  system  expressed  by 
equations  (5)  and  (6)  with  the  initial  state  2£(tQ)  and  the  input  u(t)  for 
t > t , we  seek  the  output  y(t)for  t t0.  Since  the  input  u(t)  is  expressed 
in  the  form  of  sampled  data,  both  x(t)  and  y(t)  are  determined  only  for 
discrete  multiples  of  the  simulation  sampling  interval  T. 

When  u(t)  is  a piecewise- linear  function  and  A is  time- invariant  wse 
obtain  an  approximate  solution  [13]  for  x(t)  as 


x(kT)  = T x [ (k  - 1)T]  + T B ^ u[(k  - 1)T]  + B * u(kT) 

£ C. 


(7) 


where  k is  an  integer.  For  the  A matrix  of  equation  (6),  the  matrix 
exponential  eAT  is  expressed  in  closed  form  as  follows  [14]: 
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A T 
e— 


(8) 


-aT 

e 


a 


0 
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Assuming  the  dead  time  d to  be  equal  to  mT,  where  m is  an  integer,  we  can 
express  the  input  function  u in  terms  of  the  valve  output  u^. 


u(kT)  = uv[(k  - m)T]  - u^[(k  - m - n)T] 

The  state  equation  for  simulating  the  linear  process  is  thus: 


(9) 


xx(kT) 

_ , -aT 

-aT  1 - e 

e 

a 

xx[(k  - 1)T] 

X ,(kT) 

0 1 

x2[(k  - 1)T] 

- _ 

-aT 

. -aT " 

1 - e 

- - 

a 

0 

u[(k  - 1)T]  + 

0 

0 

1 

KT 

KT 

2 

2 

L 

- 

u(kT) 


(10) 


The  process  output  y is  given  by 


y(kT)  = x1(kT) 


(ID 


Equations  (9)  through  (11)  are  used  for  the  linear  process  simulation  with  the 
piecewise  linear  input  u which  is,  in  turn,  the  change  of  valve  opening.  For 
given  T,  the  parameters  describing  the  process — a,  K,  and  m — are  time-invariant. 
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2.2  Valve  Simulation 


The  valve  works  as  a final  controL  element.  For  the  simulation  of  the  valve 
operation,  we  assume  a non-linear  valve  opening  with  respect  to  the  controller 
output  signal.  The  response  of  the  valve  is  shown  in  figure  4.  A possible 
hysterisis  of  the  valve  action  is  ignored. 

2 . 3 Controller  Simulation 

As  a simple  approach,  we  take  a Pi-controller . An  ideal  Pi-controller  [15]  has 
a transfer  function  for  a continuous  system  represented  by 


U(s) 

E(s) 


= K (1  + 
P 


(12) 


where  U(s)  and  E(s)  are  Laplace  transforms  of  the  controller  output  signal  uc, 
and  the  error  signal  e,  respectively.  Kp  and  Tj  are  the  proportional  gain 
of  the  controller  and  the  integral  time  constant,  respectively. 


Because  we  deal  with  sampled  data,  we  perform  the  z-transformation  on  equation 
(12)  to  obtain  in  the  first  order  approximation 


U(z)  = K T'z 

E(z)  p 11  z - 1 


(13) 


or 


K T ' 

U(z)  - z"1  U(z)  = K E(z)  - z_1K  E(z)  + E(z)  (14) 

P P Ti 

where  T’  is  the  controller  sampling  interval. 

The  difference  equation  corresponding  to  equation  (14)  is  thus 
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(15) 


T* 

u (kT")  = u [(k  - 1) T ] + K (1  + — )e(kT)  - K e[(k  - 1)T'] 
c c p T p 


Equation  (15)  represents  an  algorithm  for  the  discrete  Pi-controller  with  the 
parameters,  and  T^,  to  be  determined  a priori. 

The  choice  of  appropriate  values  of  and  is  commonly  [16]  made  in 
accordance  with  the  Ziegler-Nichols  guidelines  [17],  For  a single  input  single 
output  system  in  which  the  process  is  as  assumed  in  equation  (1),  the  Ziegler- 
Nichols  criteria  recommend: 

C 

K = 2—  and  Tt  = CTd,  (16) 

P Kd  I I 

where 

C =0.9,  and  C_  = 3.3. 

P 1 

The  values  of  Kp  and  Tj  are  constant  for  a non-adaptive  control  system.  For 
an  adaptive  controller,  Kp  and  are  updated  through  estimates  of  K and  d, 
as  explained  in  the  following  development. 


The  z-transform  corresponding  to  equation  (3)  [18]  is: 


H(z) 


Y(z) 

U(z) 


_£  — 1 
b^z  + b^z 

1 + a^z  ^ 


(17) 


where  Z is  the  next  larger  integer  than 
taking  the  integer  part). 
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(The  bracket  notation  indicates 


Thus  d = ZT'  - <$T'  for  some  6 such  that  0 < 6 < 1. 


The  coefficients  in  equation  (17)  are 


a 


1 


-aT 

-e 


(18) 


b,  - -a  - 
1 a 


-a6T' 
e ) 


(19) 


b 


2 


K,  -a<ST' 
— (e 
a 


-aT' 

e 


) 


(20) 


In  the  time  domain,  the  process  output  can  be  represented  from  equation  (17)  as 


y(kT')  = -ax  y [ (k  - 1)T']  + b^Kk  - A)T']  + b2u[(k  - Z - 1)T'].  (21) 


We  can  thus  write 


y(kT')  = 0T  0 


(22) 


where 


0.  = [9-,,  92’  93^T  = ^ai»  bi»  b2^ 


(23) 


0 = I0r  02,  031T 

= { -y [ ( k - l)T'l,  u[(k  - Z)T'],  u[(k  - l - 1)T']}  T (24) 
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Since  the  parameters,  a^,  b^,  and  b^,  are  to  be  determined  adaptively,  a,  6 , 
and  K are  time-dependent.  The  expressions  of  a,  6,  and  K in  terms  of  0^,  0^, 
and  0^  are  obtained  from  equations  (18),  (19),  (20),  and  (23). 


a = - 


to(-01) 


K = 


a(©2  + 0^) 

1 + 0. 


6 Jln(^  - 0 ) 

aT  K 1 


(25) 


(26) 


(27) 


Using  the  integer,  £,  which  is  given  by 


we  write  equations  (16)  in  terms  of  l and  6 as: 


C 

p K 1 " T 

Ti  - cx  {[I]  + 1 - «}  T' 

The  controller  outputs  are  determined  by  equation  (15)  with  equations  (25) 
through  (29)  when  the  parameters  0^,  0^,  and  63  are  obtained  by  the  estimation 
routine. 


(28) 


(29) 


2 . 4 Parameter  Estimator 

To  obtain  the  needed  estimates  of  these  system  parameters,  we  use  the  well-known 
method  of  recursive  least  squares  estimation  [5,  6,  9].  The  recursive  least 
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squares  algorithm,  under  quite  general  conditions,  minimizes  the  loss  function 
given  by: 


0 < A < 1 


(30) 


where  £ is  the  residual  between  the  estimates  of  the  system  parameters  of  time 
j and  their  true  values.  The  use  of  A < 1 gives  lower  weight  to  less 
recent  data.  To  implement  the  recursive  least  squares  algorithm,  we 
need  to  use  a different  form  of  process  model,  namely  the  auto-regressive 
moving  average  (ARMA)  time  series  model: 


y(t)  + a,y(t  - 1)  + . . . + a y(t  - m)  35  b.u(t  - d)  + b«u(t  - d - 1)  + . . 
i m l z 


where  y(t)  is  the  process  output  at  time  t and  u(t)  is  the  process  input  at 
time  t. 

Note  that  m and  n in  equation  (31)  are  integers,  and  the  transport  delay  (dead 
time)  is  d.  In  more  compact  notation,  a regressor  vector  and  a parameter 
vector  defined  as: 


+ b u(t  - d - n + 1) 
n 


(31) 


e = re  e 

- 1 1*  2 


e 


T 


(32) 


m+n 


T 


0 = [0r  02 


0 


T 


m+n 
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= [ — y ( t - 1),  -y(t  - 2),  , -y(t  - m)  , u(t  - d),  u(t  - d - 1) , (33) 

. . . , u(t  - d - n + 1) ]T, 

allow  the  model  given  by  equation  (31)  to  be  written  as  (Recall  equations 
(23)  to  (27)): 

y (t)  = ^T0.  (34) 

The  recursive  least  squares  parameter  estimation  algorithm  [9]  is  given  by 

0(t  + 1)  = 0(t)  + P(t  + l)0(t  + l)e(t  + 1)  (35) 


where 


£(t  + 1)  = y(t  + 1)  - 0T(t)0(t  + 1)  (36) 

P(t  + 1)  = [P(t)  - P(t)0(t)R(t)0T(t)P(t)]/X  (37) 

R(t)  = [X  + 0T(t)P(t)0(t)]_1.  (38) 

If  X is  equal  to  one,  the  magnitude  of  the  trace  of  the  P^  matrix  decreases 
monotonically.  If  X is  less  than  one,  the  trace  of  increases.  The  gain  of 
the  estimator  thus  depends  partially  upon  X,  which  is  determined  empirically. 

The  properties  of  the  _P  covariance  matrix  as  a function  of  X are  discussed 
extensively  by  Astrom  [9].  He  points  out  possible  instability  of  the  P matrix 
and  suggests  a number  of  remedies.  One  of  his  suggestions  is  adapted  in  our 
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implementation  of  the  recursive  least  squares  algorithm.  When  the  value  of 
EP_(t)j9(t)  is  less  than  a given  small  value,  the  updating  of  the  _P(t)  matrix 
is  stopped. 


To  obtain  a system  model  in  the  required  form,  so  that  the  vectors  0^  and  0^ 
are  properly  identified,  we  must  return  to  equation  (3).  However,  the  esti- 
mation algorithm  observes  the  system  inputs  and  outputs  only  at  the  controller 
sampling  interval  T" . Thus,  the  continuous  time  system  to  be  modelled  has 
Laplace  transform: 


H(s) 


_T's 

(X  - e S)(Ke 
s(s  4-  a) 


(39) 


3 . COMPUTER  SIMULATION 


3.1  Procedure 

Using  the  mathematical  model  of  the  process,  we  performed  computer  simulation 
studies,  combining  both  a constant  Pi-controller  and  an  adaptive  Pi-controller . 
In  our  simulations  of  the  adaptive  controller  we  assumed  initial  conditions 
for  the  covariance  matrix  P and  the  parameter  vector  _0  as 


P(0) 


10  0 
0 10 
0 0 1 


and  0(0)  = [-1,0, 0]T. 


(40 


Steady  state  values  were  used  as  initial  conditions  of  the  process. 


We  employed  two  different  forgetting  factors  X;  one  which  was  less  than  unity 
to  speed  up  the  convergence  and  the  other  factor  which  was  unity  to 
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keep  the  covariance  matrix  from  instability.  After  beginning  the  simulation 
or  setting  a new  reference  value,  we  used  X ■ 0.98  for  30  sampling  periods  for 

fast  convergence  of  the  parameter  estimates.  However,  since  the  trace  of  the 

covariance  matrix  grew  rapidly  when  X was  equal  to  0.98,  we  switched  to  X = 1.0 
subsequently,  which  served  to  stabilize  the  covariance  matrix. 

Although  it  is  desirable  to  take  one  X-value  which  gives  stable  values  of  I>,  we 
reinitialized  the  simulation  with  new  sets  of  initial  conditions  every  24  hours 

As  seen  in  equations  (25)  and  (27),  _0  is  one  of  the  arguments  of  a logarithm 
function.  To  ensure  positive  values  of  a and  K,  and  to  satisfy  the  boundedness 
of  6 we  imposed  limits  on  the  parameter  estimates  to  assure  that: 

<o.  e2  > o,  e3  > o. 

We  also  examined  the  values  of  £ P,  ,(t)0  (t).  If  the  growth  rate  of  the 

summed  quantity  of  P^t)?)(t)  was  less  than  and  equal  to  1.0  percent,  we  stopped 
the  parameter  update.  When  the  growth  rate  became  more  than  1.0  percent,  we 
resumed  the  parameter  update. 

To  initialize  the  controller,  we  needed  to  have  estimates  of  the  dead  time  d, 
the  process  time  constant  and  the  process  gain  K.  The  valve  response 
characteristics  were  determined  prior  to  simulation  as  shown  in  figure  4. 

As  noted  above,  two  sampling  intervals  were  required,  one  for  the  process 
and  the  other  for  the  controller  and  estimator.  We  took  one  minute  interval 
for  controller  and  estimator  sampling  and  one  second  interval  for  process 
sampling.  The  controller  sampling  time  should  be  equal  to  or  less  than  the 
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system  dead  time  for  good  stability.  (Recall  the  assumption  on  d for  equa- 
tion (17).) 

At  first,  we  simulated  the  process  with  a constant  outdoor  air  temperature  of 
-9.7°C  (14.6°F)  for  given  process  time  constant,  gain,  and  dead  time,  and 
assumed  that  there  was  no  return  air.  A step  change  in  setpoint  was  the  only 
disturbance.  From  this  simulation,  we  observed  the  behavior  of  the  parameter 
estimates  and  the  covariance  matrix  and  the  influence  of  various  values  for 
the  forgetting  factor  X. 

Next  we  considered  the  outdoor  air  temperature,  as  varying  minute  by  minute. 

We  took  a typical  daily  cycle  of  mean  outdoor  air  temperature  TQA  [19]. 

The  maximum  and  minimum  air  temperatures  were  taken  as  0°C  (32°F)  and  -11.1°C 
(12°F)  respectively. 

We  considered  the  impact  of  a stochastic  component  in  the  disturbance  measurement 
by  adding  in  pseudo  random  numbers  v.  A general  form  of  the  disturbance  w can  be 
given  simply  as: 

w = (l  - 3)tr  + 3tqa  + v (41) 

where  3,  TR  and  T^A  are  the  ratio  of  outdoor  air  to  supply  air,  the  return  air, 

and  outdoor  air  temperatures,  respectively.  In  our  simulations,  we  fixed  3 to 

be  1.0,  which  means  that  100  percent  outdoor  air  was  used.  The  resultant  process 
output  was  thus  the  sum  of  the  process  output  due  to  controller  action  and  that 
due  to  the  disturbance  w.  The  difference  of  the  resultant  process  output  from 

the  set  point  value  was  the  error  signal  that  drove  the  controller. 
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The  PI  controller  output  u£  was  initially  set  to  0.5  and  bounded  in  the 
interval  of  [0,1].  The  controller  output  u£  is  the  input  to  the  valve.  The 
variation  of  the  valve  output  is  the  input  to  the  linear  process,  while  the 
controller  output  is  the  input  of  the  estimator. 


3 . 2 Results  and  Discussion 

The  purpose  of  the  computer  simulations  was  to  test  the  validity  of  the  process 
modeling,  to  compare  the  performance  of  the  adaptive  control  algorithm  with  the 
constant  Pl-algor ithm,  and  to  provide  a basis  for  implementation  of  the  algorithm 
in  hardware  for  direct  digital  control.  We  observed  the  system  responses  due  to 
setpoint  change,  outside  air  variations,  and  change  of  process  parameters. 


Figure  5 shows  process  outputs  both  with  adaptive  Kp  and  Tj-  and  with  non- 
adaptive  and  T^..  The  initially  assumed  process  here  was 


gh(s) 


35  e 


-1.0s 


s + 0.5 


(42) 


The  parameters,  and  T^  of  the  constant  Pi-controller , were  obtained  using 
the  actual  values  of  K,  a,  and  d of  Lhe  process.  These  were  = 0.0147  and 
Tj  = 3.3.  The  adaptive  controller  needed  only  one  value,  the  dead  time, 
prior  to  simulation.  We  can  use  initial  values  of  0_  and  P^  for  all  cases  if 
plant  characteristics  do  not  change  drastically. 

We  presume  that  the  dead  time  can  be  estimated  adequately  in  advance  and  does 
not  change  drastically  under  operating  conditions.  Mean  outdoor  air  tempera- 
ture varied  during  the  simulation  period  from  0:00  a.m.  to  8:00  a.m.,  and 
a stochastic  component  with  a standard  deviation  of  0.14°C  (0.25°F)  was 
added  to  the  mean  temperature. 
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As  seen  in  figure  5,  the  adaptive  algorithm  took  about  20  minutes  to  settle 
down  due  to  an  initial  learning  period.  As  time  went  on,  the  parameter 
estimates  tended  to  the  actual  values.  The  constant  Pi-controller  damped 
more  quickly  than  the  adaptive  controller  initially,  but  after  120  minutes, 
the  process  output  of  the  adaptive  controller  had  approached  the  reference 
value  of  14.4°C  (58°F)  while  the  output  of  the  non-adaptive  controller  was 
still  oscillatory. 

The  set  point  was  changed  to  17.8°C  (64°F)  at  180  minutes  after  start-up.  Both 
outputs  converged  quite  well  to  the  new  set  point.  At  the  300  minute  mark, 
the  process  parameters  were  changed;  a from  0.5  to  0.6  and  K from  35.0  to  40.0. 
This  time  when  the  set  point  was  changed  from  17.8°C  (64°F)  to  14.4°C  (58°F), 
both  responses  were  less  stable.  However,  even  though  the  output  of  the  con- 
stant Pi-controller  remained  oscillatory  for  a 120  minute  period,  the  output 
of  the  adaptive  controller  stabilized  within  40  minutes.  From  visual  inspec- 
tion one  can  see  that  the  adaptive  controller,  in  general,  is  more  effective 
than  the  non-adaptive  controller.  We  note  that  a discrete  change  in  process 
parameters  is  quite  artificial.  It  is  done  here  merely  to  illustrate  that  the 
constant  controller  is  very  dependent  on  being  correctly  tuned. 

When  the  valve  operation  was  in  a linear  region  of  its  response  curve, 
simulations  were  reasonably  good.  But  we  obtained  unpredictable  results  when  the 
valve  response  was  nonlinear.  We  observed  nonlinear  operation  during  a start- 
up period  and  a set  point  change  period  because  of  large  disturbances.  When 
this  learning  period  was  over,  the  simulation  was  in  a linear  region.  For  the 
constant  Pi-controller , we  assumed  a linear  valve. 


17 


In  figure  6,  we  compare  the  process  outputs  controlled  by  the  adaptive  controller 

and  those  controlled  by  the  non-adaptive  controller  for  process  dead  time  of 

3.0  minutes.  The  constant  Pi-controller  was  assumed  to  be  badly  tuned  with 

Kp  = 0.0009  and  = 9.9,  which  were  computed  based  upon  a = 0.01  and  K = 200. 

Although  these  estimated  values  were  quite  extremet,  it  was  interesting  to  see 

how  much  off-set  would  occur  under  such  extreme  cases.  All  other  conditions 

were  identical  with  those  in  figure  5,  except  for  the  initial  values  of  Kp  and 

Tj.  The  adaptive  algorithm  started  with  the  initial  condition  given  by 

equation  (40).  The  Kp  and  Tj  of  the  adaptive  Pi-controller  converged  to  0.0086 

and  11.75,  respectively,  at  the  180  minute  mark,  while  those  of  the  non-adaptive 

controller  remained  as  K = 0.0009  and  Tt  = 9.9. 

P I 

Until  the  360  minute  point  in  figure  6,  the  outputs  of  the  constant  Pi-controller 
had  considerable  off-set.  But  after  360  minutes,  where  the  set  point  was 
changed  from  17.8°C  (65°F)  to  14.A°C  (58°F) , the  off-set  became  quite  small 
while  the  outputs  of  the  adaptive  algorithm  also  had  noticeable  errors  when  the 
process  parameters  were  varied.  The  adaptive  algorithm  converges  slowly  when 
a change  in  the  process  characteristics  takes  place. 

Figures  5 and  6 also  illustrate  the  performance  of  the  adaptive  control  under 
changes  in  set  point.  The  experiments  showing  performance  under  dramatic  changes 
in  process  characteristics  show  the  robustness  of  the  adaptive  algorithm:  even 

though  such  parameter  changes  can  only  happen  very  gradually  (e.g.,  seasonally) 
the  adaptive  controller  can  handle  them  quite  effectively. 
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4. 


DIRECT  DIGITAL  CONTROL  EXPERIMENTS 


4 . i Procedure 

Using  the  adaptive  controller  described  above  for  direct  digital  control  we 
performed  experiments  on  an  air  handling  unit  in  a general  purporse  laboratory  build- 
ing at  the  National  Bureau  of  Standards  (NBS)  , Gaithersburg,  Maryland.  The  control 
of  the  air  handling  unit  differed  somewhat  from  the  computer  simulation  model. 

The  valves  of  the  air  handling  unit  were  actuated  pneumatically.  In  the  simula- 
tion model,  the  valve  opening  for  the  heating  coil  was  governed  directly  by 
the  control  signal  as  if  the  valve  was  operated  depending  directly  upon  the 
electric  voltage  level.  The  apparatus  has  been  described  elsewhere  [20,21] 
including  the  microprocessor-based  controller,  the  digital-to-pneumatic  interface, 
and  the  associated  software. 

The  test  system  with  adaptive  control  is  depicted  in  figure  7.  Internal  to  the 
controller,  the  signal  u^  is  produced  based  upon  the  error  e and  the  estimated 
parameters  0_.  It  is  stored  for  d sampling  units  to  produce  the  controller 
output 

Au  = u (t  - d)  - u (t  - d - 1) . 
c c c 

The  quantity  Au  is  multiplied  by  a scale  factor  S to  become  the  operating  time 
c r 

At  of  the  motorized  actuator  during  a sampling  period.  The  actuator  regulates 
a 

the  branch  pressure  p in  the  pneumatic  system,  which  modulates  two  steam  valves 
for  the  pre-heat  coils,  one  chilled  water  valve,  and  the  dampers  for  the 
inlet,  return,  and  exhaust  air.  To  avoid  complexity,  our  experiments  were 
made  when  the  damper  opening  was  in  a steady-state  condition.  Hence  dampers 
are  not  shown  in  figure  7. 
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The  difference  in  branch  pressure  Ap  causes  a valve  opening  change  of  Au  , which 
could  not  be  monitored.  The  actual  valve  opening  u^  determines  the  flow  rates  of 
the  hot  steam,  chilled  water,  or  both.  The  air  passing  through  the  coils  com- 
bines with  the  disturbances  w to  become  the  supply  air  with  temperature  y.  The 
supply  air  temperature  is  measured  at  the  inlet  of  the  supply  air  fan,  which 
operates  continuously. 


The  motorized  pressure  regulator  (actuator),  valves,  and  coils  comprise  a 
process  described  in  equation  (1) . Thus  the  operating  time  of  the  actuator 
is  an  input  to  the  process.  Nonlinearity,  changes  in  gain  and  time  constant, 
changes  in  dead  time,  and  hysteresis  in  the  process  were  all  observed.  An 
effective  adaptive  controller  can  be  expected  to  compensate  for  many  different 
changes  in  the  process,  provided  they  do  not  occur  too  rapidly. 

With  experience,  it  was  possible  to  improve  the  effectiveness  of  the  control 
algorithm.  Unlike  the  computer  simulation,  the  direct  digital  control  on  the 
real  hardware  had  a narrow  range  of  operation  which  did  not  allow  the  supply 
temperature  to  be  raised  by  more  than  24°C  (75.2°F),  which  had  been  set  as  a 
high  limit. 

Thus,  restrictions  on  the  controller  output  and  the  parameter  estimates  were 

needed  such  that  u € -u  ,u  1 and  0 €(-1,0).  The  limiting  value  u was  computed 

C I R R ] i K. 

from  the  following  expression: 


u 


R 


K (1  + 
p,max 


) e 

1,0 


(43) 


where 
and  T 


'R 


1,0 


is  the  operating  temperature  range,  is  the  maximum  value  of  K^, 


is 


the  integral  time  given  by  equation  (29)  with  5=0.  Based  upon 
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previous  work  [20],  we  set  K qv  to  0.1,  and  e to  10.0°C.  Equation  (A3)  was 

‘ PlIuaX  ^ 

taken  to  have  a similar  form  of  equation  0-5),  but  it  was  a choice  of  convenience. 

Furthermore,  the  coefficient  Cp  in  equation  (28)  and  in  equation  (29) 

were  set  to  0.5  K_.,  and  1.1,  respectively. 

P max 

The  scale  factor  Sp  is  an  empirical  conversion  of  the  change  in  control  signal 
into  the  actuator  operating  time.  It  proved  to  be  of  great  importance  in  our 
tests.  The  following  relationship  was  used: 

Ata  = -SfAuc  + 4 (44) 

The  constant  t^  represents  the  leftover  time  in  a sampling  period  after  the 
actuator  has  operated.  This  was  needed  since  the  motorized  actuator  could  not 
operate  for  a time  interval  of  less  than  0.1  seconds.  The  FID  operating  system 
[21]  was  capable  of  accumulating  this  leftover  time  until  it  was  greater  than 
or  equal  to  0.1  seconds.  The  constant  t^  may  be  omitted  from  the  r.h.s.  of 
equation  (AA)  for  lower  values  of  |tp/Ata|. 

In  the  parameter  estimator,  the  initial  conditions  of  the  parameter  vector  9_ 
were  set  to  0(0)  = [-0.9,  0,  0]  . The  forgetting  factor  X was  employed  in  the 
same  manner  as  in  the  computer  simulation.  After  2A  hours  of  operation,  it 
was  reinitialized  at  the  value  0.98.  A flow  chart  showing  essential  parts  of 
the  algorithm  is  given  by  figure  8.  The  computer  program  of  the  adaptive  control 
portion  is  appended. 

4 # 2 Results  and  Discussion 

Tests  were  performed  from  mid— December  1981  until  the  end  of  February  1982. 

During  this  period,  the  outdoor  temperature  varied  from  mild  to  severely  cold. 


Most  of  our  testing  was  done  to  determine  the  system  response  to  a step  change 
of  the  set  point  between  15.0°C  (59°F)  to  18.3°C  (65°F). 

For  a set  point  of  18.3°C,  figure  9 shows  time  series  of  the  process  output  y, 
the  error  e,  the  parameter  estimates  0i,  the  integral  time  Tj,  the 
proportional  gain  K^,  the  control  signal  uc,  and  the  actuator  operating  time 
At.  The  adaptive  control  algorithm  was  initialized  at  the  beginning  of  the 

a 

run.  Table  1 describes  the  operating  parameters  for  the  controller  for  this 
run. 

Figure  9 shows  that  the  parameter  estimates  converged  after  60  samples 
(i.e.,  20  min).  The  value  of  Kp  remained  at  its  given  maximum  of  0.1, 
while  Tj  varied  actively  as  0^  changed.  The  output  y took  longer,  about 
80  minutes  to  settle  at  reference  value.  Consequently,  the  control  signal 
u£  was  large  until  y became  steady.  The  signals  of  y and  Ata  are  very 
similar . 

Figure  10  represents  a typical  response  associated  with  set  point  changes. 

We  maintained  reasonably  constant  process  outputs  prior  to  the  set  point 
change.  Since  the  scale  factor  depends  on  the  characteristics  of  the 
actuator  and  its  interface  with  the  controller,  two  different  scale  factors 
Sp  were  applied  to  improve  the  system  response  characteristics.  That  is, 
when  the  set  point  was  changed,  the  scale  factor  was  changed  simultaneously, 
using  Sp  = 30.0  for  r = 15.0°C,  and  Sp  = 20.0  for  r = 18.3°C.  The  scale 
factors  were  selected  so  as  to  provide  a consistent  level  of  excitation  in 
the  output  signals.  The  selection  was  empirical,  based  on  past  observations. 
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The  scale  factor  Is  very  important  to  the  adaptive  Pi-control,  and  must  be 
estimated  carefully  for  the  best  performance. 

As  shown  in  figure  10,  there  were  overshoots  in  y when  the  reference  was 

changed  (step  change).  As  would  be  expected,  the  rate  of  damping  and  the 

amount  of  overshoot  are  strong  functions  of  S^.  Parameter  estimation  was 

quite  rapid  after  disturbances  due  to  step  changes.  Figure  11  shows  an 

example  of  good  performance  due  to  satisfactory  choice  of  S . 

r 

When  the  direct  digital  control  experiments  are  compared  with  the  computer 
simulations,  one  sees  reasonably  good  agreement  in  output  damping  characteristics. 
When  the  process  is  stable,  as  assumed  in  the  computer  simulations,  the  adap- 
tive Pi-controller  becomes  a constant  Pi-controller  once  the  learning  period 
is  over. 

5.  CONCLUSION 

Our  testing  has  shown  that  a self-tuning  Pi-controller  employing  the  recursive 
least  squares  estimator  can  be  used  satisfactorily  for  control  of  an  air  handing 
unit.  When  the  unit  is  operating  in  a linear  region,  performance  is  particu- 
larly accurate.  Even  in  nonlinear  operation  regions,  performance  is  generally 
better  than  that  of  a constant  parameter  controller.  Moreover,  although  the 
direct  digital  controller  that  was  implemented  had  a different  interface 
between  the  valve  and  controller  than  the  mathematical  model,  good  agreement 
was  observed  between  experiment  and  simulation. 
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For  the  purposes  of  commercial  implementation,  the  most  important  result  of 

our  testing  was  the  observation  of  the  importance  of  the  sacle  factor,  S , 

r 

which  may  be  interpreted  as  the  product  of  the  component  gains  in  the  system. 
It  was  seen  that  the  correct  choice  of  the  scale  factor  permits  greater  effi- 
ciency of  operation  of  the  controller  and  parameter  estimator.  More  research 
is  needed  to  develop  a procedure  that  includes  a method  for  changing  the  scale 
factor  because  the  gain,  time  constant,  and  dead  time  of  HVAC  systems  depend 
on  the  supply  air  temperature.  These  dependencies,  along  with  non-linearities 
in  the  system  are  not  included  in  the  system  model  discussed  in  this  report. 
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Table  1.  Input  Data  for  figure  9 


Item 

Symbol 

Value 

Unit 

set  point 

r 

18.3 

°C 

sampling  period 

T" 

20.0 

sec 

process  dead  time 

d 

21.0 

sec 

operating  temperature  range 

eR 

10.0 

°C 

scale  factor 

SF 

20.0 

- 
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Figure  1.  A simple  heating,  ventilation,  and  air  conditioning 
(HVAC)  system  - air  handling  unit 


Figure  2.  Schematic  diagram  of  the  self-tuning  controller  system 
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Figure  3. 


Block  diagram  of  the  assumed  model  of  the  process 
to  be  controlled 


Figure  4.  The  valve  response  curve 
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Figure  5-*  Simulated  process  output  assuming  a dead  time  of  1.0  minute; 

solid  and  dotted  lines  represent  adaptive  and  nonadaptive 
control,  respectively 
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Figure  6.  Simulated  process  output  assuming  a dead  time  of  3.0  minutes; 

solid  and  dotted  lines  represent  adaptive  and  nonadaptive 
control,  respectively 


31 


Figure  7.  Schematic  diagram  of  the  direct  digital  control  system 
with  the  adaptive  control  algorithm 
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Figure  8 
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A flow  chart  of  the  adaptive  algorithm  for  direct  digital  control 
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Figure  9.  Experimental  results  of  direct  digital  control  at  the 
set  point  18.3°C 
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Figure  10.  Experimental  results  of  direct  digital  control  with  two 
set  points  15.0°C  and  18.3°C  (typical  case) 
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Figure  11,  An  example  of  good  performance  due  to  satisfactory 
choice  of  Sp 
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APPENDIX 
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*******************************************:************************* 
PICNTA:  THE  ADAPTIVE  P I -CONTROLLER 

VERSION  0628A  C.P. 

THIS  PICNTA  SUBROUTINE  IS  CALLED  BY  THE  MAIN  PROGRAM  FID09B.FOR 
AND  CALLS  THE  ADAPTIVE  PI -CONTROL  ALGORITHM. 

THE  FOLLOWING  DATA  ARE  ENTERED  VIA  FID09B  PROGRAM: 

ER=  TEMPERATURE  RANGE  IN  C 

TDED=  TRANSPORTATION  DELAY  TIME  (DEAD  TIME)  IN  S 

SCALE2  SCALE  FACTOR  ( SYSTEM  GAIN) 

DTC=  SAMPLING  PERIOD  IN  S 

DUMAX2  MAXIMUM  ACTUATOR  OPERATION  TIME  IN  S 
RIN2  INPUT  DEVICE 

ROUT2  OUTPUT  DEVICE 

HYSTIM2  COMPENSATION  TIME  FOR  HYSTERISIS  IN  S 
TREF=  SET  POINT  TEMPERATURE  IN  C 
ANAI=  SUPPLY  AIR  TEMPERATURE  IN  C 
REST2  LEFTOVER  TIME  AFTER  THE  ACTUATOR  RESPONDED 


SUBROUTINE  P ICNTA( N, DU, EO, REST) 

INTEGERS  N 

INTEGER  DUMAX, RIN, ROUT 
REAL  KP 
REAL*8  ANA I 

COMMON  /SETPNT/  TREFC  16) 

COMMON  /REG/  ER(  16) , TDED(  16) ,SCALE(  16) ,DTC(  16) ,DUMAX(  16) , 
8 RIN( 16) ,ROUT( 16) 

COMMON  /ANALOG/  ANAK16) 

COMMON  /BK3/  THETA(  3) , PHI(  3) , P( 3, 3) 

G 

C YMAX  2 HIGHEST  PROCESS  OUTPUT  VALUE  IN  C 

C 

DATA  YMAX, IH IT/22. 0,0/ 

C 

CALL  DISINT 
J=RIN(  N) 

C 

TCNT=DTC(N)/60. 

REF=TREF( N) 

YOUT= ANAI(  J) 

E0=  REF- YOUT 
ERANGE2 ER( N) 

TDEAD=TDED(N)/60. 

TACT=SCALE(  N) 

DUMAX 1 2 DUMAX(  N) 

C TO  PREVENT  FROM  REACHING  HI-CUTOFF  POINT 

IF< YOUT. GE.YMAX+ 1.0)  GOTO  10 
IF(  YOUT. LE. YMAX)  GOTO  20 
IF( IHIT.GT.2)  GOTO  20 
10  DU=DUMAX1 

I HIT2 IHIT+1 
WRITE( 5 , 100) 

100  FORMAT!  2X, ’ ***  TOO  HIGH  TEMP  ********’) 

GOTO  30 
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1580 

1590 

1600 

1610 

1620 

1630 

1640 

1650 

1660 

1670 

1680 

1690 

1700 

1710 

1720 

1730 

1740 

1750 

1760 

1770 

1780 

1790 

1800 

1810 

1820 

1830 

1840 

1850 

1660 

1870 

1880 

1890 

1900 

1910 

1920 

1930 

1940 

1950 

I960 

1970 

1980 

1990 

2000 

2010 

2020 

2030 

2040 

2050 

2060 

2070 

2080 

2090 

2100 

2110 

2120 

2130 

2140 

2150 


20  IHIT=0 

CALL  AC AP I ( TCNT , TDEAD , REF , YOUT , UCNT , DELU , I COURT , ERAN GE , T I C , KP ) 
C 

C ACTUATOR  OPERATIRG  TIME 

C 

DU=  - DELU*TACT+ REST 


C 

WRITE! 5 , 200)  I COUNT, (THETA! I) , 1=1,3) , KP, TIC, UCRT, E0 , DU,  YOUT 
200  F0RT1AT(  15 , 3F8. 4 , 3F 10 . 4 , 3F7 . 3) 

30  CALL  ENA I NT 

C 

RETURN 

END 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


************************************************************************** 
ACAPI:  ADAPTIVE  CONTROL  ALGORITHM  FOR  PI-CONTROLLER 
TCNT=  SAMPLING  PERIOD  IN  MIN 

TDED=  DEAD  TIME  IN  MIN  — TDED  SHOULD  BE  GREATER  THAN  TCNT 
REF=  SET  POINT  (REFERENCE  VALUE) 

YOUT=  PROCESS  OUTPUT 

UCNT=  BOUNDED  CONTROLLER  OUTPUT 

DELU=  INPUT  TO  PROCESS 

I COUNT5  COUNT  INDEX 

ERANGE=  PROCESS  OPERATING  RANGE 

TIC=  INTEGRAL  TIME  IN  MIN 

KP=  PROPORTIONAL  GAIN 

*************************************************************************** 


SUBROUTINE  AC AP I ( TCNT, TDED, REF, YOUT, UCRT, DELU, ICOUNT, ERANGE, 
8 TIC , KP) 

REAL  LAMBDA, LAM(  2) ,KP 
DIMENSION  UTEMP( 20) 

COMMON  /BK1/  UCN(  20) .ERROR, ERL AST, NDEDTC , URANGE 
8 /BK2/  YEST, LAMEDA, PPHSUM, PPHSM1 , ILSET, NLSET, EPSPH 

8 /BK3/  THETA(3) ,PHI(3) ,P(3,3) 

DATA  EPSREF, EPSPH, LAM(  1) ,LAM( 2) /0. 01 ,0.01 ,0.98, 1.00/ 

C ITMAX=24  HRS*60/( 20/60) 

DATA  ITCNT, ITINTL, I TMAX/0 , 30 , 4320/ 

C 

C INITIALIZATION 

C 


IF(  ITCNT. GE. 1)  GOTO  30 
DO  10  1=1,3 
DO  10  J= 1 ,3 
THETA(  I) =0.0 
P(  I, J)=0.0 

I F ( I.EQ.J)  P(  I,J)  = 1.0 
10  CONTINUE 


THETA( 1 ) =-0.9 
C 

NDEDTC= I F I X( TDED/TCNT) 

NP=NDEDTC+5 

REF  1 = REF 

NLSET=  30 

LAMBDA5 LAM(  1) 

ERL  AS  T=  0.0 
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2160 
2170 
2180 
2190 
2200 
2210 
2220 
2230  20 
2240 
2250 
2260  C 
2270  C 
2280  G 
2290  30 
2300 
2310 
2320  C 
2330  C 
2340  C 
2350  40 
2360 
2370  50 
2380 
2390  60 
2400 
2410  G 
2420  C 
2430  C 
2440 
2450 
2460 
2470 
2480  C 
2490  C 
2500  C 
2510 
2520 
2530  70 
2540 
2550  80 
2560  C 
2570  C 
2580  G 
2590 
2600 
2610  C 
2620  C 
2630  C 
2640 
2650 
2660 
2670 
2680 
2690  C 
2700  C 
2710  C 
2720  90 
2730 


PPHSM12 1.0 

IF(  ABS(  YOUT-REF) . GT. 0 . 02*ERANGE)  YOUT= 1 . 01*REF 

ERROR2 REF-YOUT 

URANGE=0. 1 

DELU=0 . 0 

YEST1=0. 0 

DO  20  J= 1 , NP 

UCN( J) =0.00 1 

I LSET=  0 

GOTO  90 

STEP  CHANGE  IN  SET  POINT 

IF(  ABS( ( REF-REF 1 ) /REF 1 ) . LE.EPSREF)  GOTO  40 

ILSET=0 

LAMBDA2  LAM(  1) 

CONTROLLER  OPERATION 

ERROR2 REF-YOUT 

I F ( ITCNT.GT. 1TINTL)  GOTO  60 

ERROR2 -ERROR 

YOUT2 Y0UT-ERR0R*2 . 0 

CALL  P IZN( TCNT, UCNT, ERANGE, TIC, KP) 

UGNT0=UCNT 

BOUNDING  THE  CONTROLLER  OUTPUTS 

ULO=-URANGE 
UE I 2 URANGE 

IF(UCNT.LT.ULO)  UCNT=ULO 
IF(UCNT.GT.UHI)  UCNT=UHI 

SHIFTING  ONE  TIME  STEP  BACKWARDS 

UTEMP( 1) 2 UCNT 
DO  70  J 2 2 , NP 
UTEMP ( J ) 2 UCN ( J- 1 ) 

DO  80  J2 1 , NP 
UCN(  J) = UTEMP ( J) 

INPUT  SIGNAL  TO  THE  PROCESS 

IF(  I TCNT. LE. ITINTL)  GOTO  90 
DELU2 UCN( NDEDTC+ 1 ) -UCN( NDEDTC+2) 

DEALING  WITH  SATURATION  CONDITION 

I F ( UCNT. EQ. ULO. AND. UCNT. EQ. UCN( 2) ) DELU2UCNT0-UCNT1 

I F ( UC  NT . EQ . ULO . AND . DELU . GT . 0 . 0 ) DELU2  0 . 0 

IF(  UCNT. EQ. UHI . AND. UCNT. EQ. UCN( 2) ) DELU=UCNT0-UCNT1 

IF(  UCNT. EQ. UHI . AND. DELU. LT. 0 . 0)  DELU=0.0 

IF( UCNT. EQ. UCN(  NP) ) DELU=5.*DELU 

ESTIMATION  OF  PARAMETERS 

YEST2 YOUT-REF 
PHI ( 1) =-YESTl 
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2740 
2750 
2760 
2770 
2780  C 
2790  C 
2800  C 
2810 
2820 
2830 
2840 
2850 
2860  C 
2870 
2880 
2890 
2900 
2910  C 
2920  C 
2930  C 
2940 
2950 
2960  C 
2970 
2980 
2990  C 
3000  C 
3010  G 
3020  C 
3030  C 
3040  G 
3050 
3060 
3070 
3080 
3090 
3100  C 
3110  C 
3120  C 
3130 
3140 
3150 
3160 
3170  C 
3180  C 
3190  C 
3200 
3210 
3220 
3230 
3240 
3250 
3260 
3270 
3280 
3290  C 
3300 
3310 


PH I ( 2 > = UCN( NDEDTC+ 1 ) 

PHI ( 3) = UCN( NDEDTC+2) 

IF(  ILSET. GT. NLSET)  LAMBDA=LAM(  2) 

CALL  ESTIM 

CHANGE  CURRENT  VALUES  TO  THE  LAST 

YEST1= YEST 
REF  1 = REF 
ERLAST= ERROR 
PPBSril  = PPHSUM 
UCNT1  = UCNT0 

ILSET= ILSET+1 
ITCNT= ITCNT+1 

I F ( 1TCNT.EQ. ITINTL)  ILSET=0 
IF(  1TCNT.LE. ITINTL)  GOTO  50 

RE- INITIALIZATION  EVERY  24  HOURS 

I COUNT= ITCNT- ITINTL 

IF( 1C0UNT.EQ. ITMAX+1)  ITCNT=0 

RETURN 

END 

*********************************************************************** 
PIZN:  ADAPTIVE  PI-  CONTROLLER 

**********#*****************v******************************************* 

SUBROUT I NE  PIZN(  TCNT, UCNT, ERANGE , T I C , KP ) 

REAL  KP , KPMIN , KPMAX 

COMMON  /BK1/  UCN(  20) .ERROR, ERLAST, NDEDTC, URANGE 
8 /BK3/  THETA( 3) , PHI( 3) , P( 3,3) 

DATA  CTI/1. 1/.KPMAX/0. 1/, KPMIN/0. 01/ 

ADAPTIVE  PI -CONTROLLER 

ALF=-ALOG( - THETA  ( 1 ) ) /TCNT 

GN= ALF*( THETA( 2) +THETA( 3) ) /<  1 . +THETA(  1)) 

DELT=-ALOG( ALF*THETA(  3) /GN-THETA(  1 ) ) /(  ALF*TCNT) 

IF(DELT.GE. 0.999)  DELT=0.999 

MODIFIED  ZIEGLER- NICHOLS  CRITERIA 

TDED 1 = ( NDEDTC+ I-DELT) *TCNT 
TIC=CTI*TDED1 
T I F I X= CT I * ( NDEDTC+ 1 ) *TCNT 
CKP=0 . 5*KPMAX 
KP=CKP/( GN*TDED1) 

1F( KP. GE. KPMAX)  KP=  KPMAX 
IF(KP.LE. KPMIN)  KP=KPMIN 

UCNT=  UCN(  1 ) +KP*(  ( 1 .+TCNT/TIC)*EIU10R- ERLAST) 

URANGE= KPMAX*  ( 1 . +TCNT/TIFIX) *ERANGE 

RETURN 

END 
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C 

C ESTIM:  ESTIMATION  OF  PARAMETERS  USING  A RECURSIVE  LEAST  SQUARES 

C METHOD 

C 

c***************************************************************************** 

c 

SUBROUTINE  ESTIM 
REAL  LAMBDA 
DIMENSION  PPHI(3) 

COMMON  /BK2/  YEST, LAMBDA, PPHSUM, PPHSM1 , ILSET, NLSET, EPSPH 
8 /BK3/  THETA(  3) , PHI ( 3) , P( 3 , 3) 

C 

DO  10  J= 1 ,3 
PPHI ( J) =0 . 

DO  10  K= 1 ,3 

10  PPHI ( J) =PPHI ( J) +P(  J, K) *PHI(  K) 

PPHSUM=0 . 

DO  20  J=  1 , 3 

20  PPHSUM=PPHSUM+PPHI(  J) 

C 

E1  = 0. 

DO  30  J= 1 ,3 

30  E1  = E1+THETA(  J)*PHI< J) 

E= YEST- El 
G 

I F( ILSET.LT. NLSET)  GOTO  40 

IF( PPHSM1 . NE. 0. . AND. ABS( ( PPHSUM-PPHSM1) /PPHSM1 ) .LE. EPSPH)  RETURN 
C 

40  CONTINUE 

DO  50  J= 1 ,3 

50  THETA( J) = E*PPHI(  J) +THETA( J) 

I F ( THETA( 1 ) . GE . 0.0)  THETAI 1) =-0.0001 
IF(THETA( 1) .LE.-1.0)  THETA(  1 ) = -0 . 999 
IF(THETA(2) .LT.O. ) THETA(  2) = 0.0001 
IF(THETA(3) .LT.O. ) THETA(  3 ) = 0.0001 
G 

R= LAMBDA 
DO  60  J= 1 ,3 

60  R=R+PHI( J)*(P( J, 1)*PHI<  1)+P(  J,2)*PHI(2)+P<  J,3)*PHI<3)) 

R=l./R 

G 

DO  70  J= 1 ,3 
DO  70  K=l,3 

70  P(J,K)  = (P(J,K) -PPHI  ( J) *R*PPHI ( K) ) /LAMBDA 

C 

RETURN 

END 
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